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Abstract. 

The characteristics of an x-ray spectrum can greatly influence imaging and related 
tasks. In practice, due to the pile-up effect of the detector, it’s difficult to directly 
measure the spectrum of a CT scanner using an energy resolved detector. An 
alternative solution is to estimate the spectrum using transmission measurements 
with a step phantom or other CT phantom. In this work, we present a new 
spectrum estimation method based on indirect transmission measurement and model 
spectra mixture approach. The estimated x-ray spectrum was expressed as weighted 
summation of a set of model spectra, which can significantly reduce the degrees of 
freedom (DOF) of the spectrum estimation problem. Next, an estimated projection 
can be calculated with the assumed spectrum. By iteratively updating the unknown 
weights, we minimized the difference between the estimated projection data and 
the raw projection data. The final spectrum was calculated with these calibrated 
weights and the model spectra. Both simulation and experimental data were used 
to evaluate the proposed method. In the simulation study, the estimated spectra 
were compared to the raw spectra which were used to generate the raw projection 
data. For the experimental study, the ground truth measurement of the raw x-ray 
spectrum was not available. Therefore, the estimated spectrum was compared against 
spectra generated using the SpekCalc software with tube configurations provided by 
the scanner manufacturer. The results show the proposed method has potential 
to accurately estimate x-ray spectra using the raw projection data. The difference 
between the mean energy of the raw spectra and the mean energy of the estimated 
spectra were smaller than 0.5 keV for both simulation and experimental data. Further 
tests show the method was robust with respect to the model spectra generator. 
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X-ray spectral characteristics greatly influence imaging and imaging related tasks such 
as Monte Carlo (MC) based radiation dose calculation, artefacts reduction techniques, 
and dual-energy material decomposition. In practice, it is difficult to directly measure 
the spectrum of a CT scanner with an energy resolved detector due to the limited 
count rate of the detector. Instead, indirect methods are widely used to estimate 
the spectrum of a CT scanner. These methods can be roughly classified into model- 
based methods which estimate the spectrum using empirical or semi-empirical physical 


models (Birch and Marshall, 1979; Tucker et al ., 1991; Boone and Seibert, 1997), 


Compton-scattering measurement (Yaffe et a/., 1976; Matscheko and Ribberfors, 1987 


Matscheko and Carlssor 

l, 1989; Matscheko et a/., 1989; Gallardo et a/., 2004; Maeda 

et al, 

2005 

), transmission measurements (Archer and Wagner, 1982 

Hassler et al, 

1998; Waggener et al, 

1999) Sidky et al, 2005; Zhang et al, 2007; ! 

Xu et al l 2009; 

Duan et al, 

2011; Lin et al, 2014) and MC simulation (Llovet et al, 

2003; Ay et al, 

2004| Mainegra-Hing and Kawrakow, 2006; Bazalova and Verhaegen, 2007; Miceli et al, 


2007a). 


The aim of this work was to develop an indirect transmission measurement-based 
spectrum estimation method using only raw projection data of a physical phantom 
with known density to estimate the corresponding spectrum of a CT scanner. To 
determine the x-ray spectrum of a CT scanner, energy-resolved detectors such as 
cadmium zinc telluride (CdZnTe) ( |Miyajima et al , 2002; Fritz et al. , 2011), cadmium 
telluride (CdTe) (Miyajiriiaj 2003| Redus et al., 2009) or high purity germanium (Fewell 


et al, [1981) are often used to directly measure the spectra. However, due to the 


limited count rate of these detectors and the high photon flux of CT scanners, it is 
not easy to directly measure the spectrum of a CT scanner because of the detector 
pile-up effect. Besides, direct measurements using these energy-resolved detectors may 
suffer from environmental requirement (such as low temperature requirement) or hole 
trapping effects, which yields low-energy tailing of the energy spectrum (Koenig et al. 


2012). Instead of direct measurement, Compton scattering techniques can be employed 


to reconstruct the spectrum ( 

Yaffe et al, 

1976; 

Matscheko and Ribberfors, 

1987; 

Bartol 


2013). In this method, an energy spectrometer is placed outside the primary ray to 


receive Compton scattering photons from a scatter object with significantly reduced 
photon events. The incident photon energy was then calculated with the corresponding 
energy of the Compton scattering photons and the scattering angle. 

Alternatively, the properties of a polychromatic spectrum can be extracted from 
projection measurements. For example, a harder spectrum will generally yield less 
attenuated projection measurements and a softer spectrum more attenuated projection 
measurements. Thus, an alternative solution for spectral measurement is to estimate the 


spectra with transmission measurements (Sidky et al, 2005; Zhang et al, 2007; Duan 


et al, 2011; Zhang and Zhang, 2013). These methods usually employ a step or a wedge 


phantom and a polychromatic forward projection equation formulated as a discrete 
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linear system. In this linear system, each energy bin of the spectrum is described 
as an unknown variable. When employing a step wedge phantom, the total number 
of unknown variables maybe larger than the number of transmission measurements. 
Thus the linear system may be ill-conditioned and the expectation-maximization (EM) 
algorithm is often used to solve the problem. In order to further improve the estimation 
accuracy of the characteristic peak of the spectrum, an improved parameter spectrum 


model was also employed to recover the spectrum (Zhang and Zhang, 2013). 


In this work, we propose a new spectrum estimation method based on indirect 
transmission measurements and model spectra mixture approach. Instead of using a 
dedicated step or wedge phantom for transmission measurements, this method aims to 
calculate the propagation path length (PPL) for each of the segmented materials for each 
detector pixel and the raw projection data (after taking a logarithm) is compared to the 
polychromatic reprojection data which is calculated using the PPL and an estimated 
spectrum. The estimated spectrum is then iteratively updated to minimize the difference 
of the raw projection data and the polychromatic reprojection data of the volumetric 
images. In order to make the iterative algorithm robust and stable, the estimated 
spectrum was expressed as the weighted summation of a set of model spectra which 
significantly reduce the degrees of freeom (DOF) of the spectrum estimation problem. 
The model spectra can be generated with either Monte Carlo simulation or other 
analytical spectrum generators. Spectra recovered with experimental phantom data 


were compared against the spectrum generated using SpekCalc software (Poludniowski 


et al, 2009) with x-ray source configurations provided by the scanner manufacturer. 


2. METHODS 


The problem of spectrum estimation based on transmission measurements is primarily 
an inverse problem. In order to determine the spectrum, we have to determine the 
content of each energy bin of the spectrum using the raw projection data. Depending 
on the width of the energy bin, the inverse problem typically has tens of unknown 
variables. Namely, there are tens of DOF of the inverse problem. Instead of estimating 
each energy bin content for a spectrum fl(E), a weighted summation of a set of model 
spectra Cli(E) is used to express Q(E). 

M 

Q(E) = Y t c i n i (E). ( 1 ) 


1=1 


Here M is the number of the model spectra and c ( is the weight on the respective model 
spectrum. The model spectra can be generated with MC simulation tool-kit or any 


other kinds of spectrum generator, such as SpekCalc (Poludniowski and Evans, 2007 


Poludniowski, 2007; Poludniowski et al, 2009). Now the number of unknown variables 


of the spectrum estimation problem can be reduced from the number of energy bins to 
the number of model spectra M, i.e. a significant reduction of the DOF of the spectrum 
estimation problem is achieved by using the model spectra mixture approach. 
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Figure 1 . A flowchart of the proposed spectrum estimation method. The method 
starts with the raw projection data. After an initial CT image reconstruction, the 
image is segmented into different components, and a polychromatic reprojection is 
performed with the segmented image and a weighted summation of a set of model 
spectra. The unknown weights are calibrated by iteratively minimizing the difference 
of the estimated reprojection data and the raw projection data. The final spectrum is 
calculated using the weights and the model spectra. 


Figure [l] shows the flowchart of the proposed indirect transmission measurements- 
based spectrum estimation method. The method starts from the original raw projection 
data. After an initial CT image reconstruction, the image is segmented to different 
components. The linear attenuation coefficients which can be obtained from the NIST 
database are assigned to these components. The polychromatic reprojection is then 
performed using Q(E) and the segmented image to generate an estimated projection 
data. By iteratively updating the weights q, we can find a set of optimal weights q to 
minimize the difference of the estimated projection data and the raw projection data. 
The final spectrum is calculated by linearly combining the model spectrum fli(E) with 
the corresponding weight q. 


2.1. Polychromatic reprojection 


The polychromatic projection of an object can be represented as 


I = N 


d E Cl(E) f]{E) exp 


— / fi(E, s)ds 


( 2 ) 


J o Jo J 

where N is the total number of photons, rj(E) is the energy dependent response of 
the detector and can be considered proportional to photon energy E because most CT 
scanners use energy-integrating detectors. E max is the maximum photon energy of the 
spectrum. fj(E, s) is the energy-dependent linear attenuation coefficient and l is the 
propagation path length for each ray. Note that I is detector pixel dependent and the 
detector channel index is omitted for convenience. With no object in the beam path, 
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the flood field Jo can be expressed as follows: 

r^max 


Io = N 


dEQ(E) r)(E). 


After applying the logarithmic operation, the projection data can be expressed as: 


p = log 


I 


= log 


jE max 


f t f max dE ft(E) r/(E) exp -J^(E,s)ds 


(3) 


(4) 


For a given spectrum estimate and segmented image, an estimated polychromatic 
reprojection dataset p can be calculated using equation Q. In addition, from the 
normalization condition, we have 

['Em ax 

/ dEQ(E) = l, (5) 


/ dEfti(E) = l. 

Jo 

Substituting equation © into equation ^ and combining with equation ([6]), 
get the following constraint condition for weights c, 

M 

5 > = 1 - 

2=1 

This constraint will be used to calibrate the weights of the model spectra. 


( 6 ) 


we can 


(7) 


2.2. Calibrating the weights 

For each detector pixel measurement, p m , the difference between the measurement 
and the corresponding estimated reprojection, p, should be minimal if the estimated 
spectrum matches the actual raw spectrum. We can minimize the difference between 
p m and p by updating the weights c. Thus the optimal set of c can be achieved by 
solving the following optimization problem, 

c = argmin \\p m - p\\\. (8) 

C 

In order to keep the solution of the above equation physically meaningful, a non¬ 
negative constraint can also be applied to the weights. Together with the constraint 
equation ([T]) , the above unconstrained optimization problem is formulated into a 
constraint problem as follows, 

M 

c = argmin \\p m - p\\l, s.t. ^ q = 1, and q > 0. (9) 

c • i 

2=1 

To solve the constraint optimization problem of equation ([9]), in this study, we have 
reformulated it into an unconstrained problem and solved the unconstrained problem 
using a simple multi-variable downhill simplex method. Thus the original constraint 
problem of equation © can be solved with a sequential optimization approach, 
minimizing the objective function, followed by normalizing the solution and enforcing 
non-negative constraint sequentially. The algorithm can be summarized as follows: 
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Algorithm (Sequential optimization) 

1. Set k=0, choose c 0 . 

2. repeat 


Ak II 2 

P \\2 


3. cr +1 = argmin \\pn 

4. c k+1 = c fc+1 /sum(c fc+1 ) 

5. c k+1 = (c fc+1 ) + 

6. k <— k + 1 

7. until stopping criterion is satisfied. 

Line 1 of the algorithm establishes the initial values of the iterative procedure. In 
this work, the weights for the initial model spectra, c 0 , were either set to be equal or to 
provide the hardest possible spectrum (i.e. the weight of the hardest model spectrum 
was set to one and all other weights were set to zeros). Lines 2 to 7 iteratively loop to 
calibrate the optimal weights beginning with the initial weights established in line 1 and 
terminating when either the maximum iteration number has been met or the difference 
between the weights of the current and previous iteration (c& and c*, + i) is smaller than e. 
In this work, the following stopping criterion was used: £ = 10 -5 . Within the loop, line 3 
is an inner update of c to minimize the quadratic error between the measured projection 
data p m and the estimated reprojection data p. This is an unconstrained optimization 
problem and a multi-variable downhill simplex method (PressJ [2007) was employed to 
solve it. The updated weights, c, are then normalized (line 4) and constrained to be 
non-negative, with weights < 0 set equal to 0 (line 5). Finally, the current values are 
assigned to be the starting point of the next iteration in line 6. The optimal weights 
generated with the above algorithm are used to determine the final spectrum. 

It should be noted that the original constraint problem equation is not 
complicated and it can also be solved directly using the methods that can intrinsically 


handle constraints, such as BFGS-B (Byrd et al, 1995) and augmented Lagrangian 


method (ALM) (Nocedal and Wright, 2006). 


2.3. Model spectrum generation 


In order to generate the model spectra, the MC simulation tool-kit Geant4 (Agostinelli 


et al ., 2003) can be employed. Geant4 provides comprehensive physics modelling 


capabilities embedded in a flexible structure and it offers a complete electron interactions 
modelling, such as the bremsstrahlung effect and characteristic radiation. With the 
low energy package G4EMLOW which can describe the electromagnetic interactions of 
photons, electrons, hadrons and ions with matter down to very low energies (eV scale), 
Geant4 can precisely model the x-ray photons generation in many kinds of targets and 


the results of energy spectra have been well validated for keV energy levels (Taschereau 


et al, 

2006; 

Miceli et al ., 

2007a[b; 


The x-ray source modelled with Geant4 includes an electron gun, a tungsten 
target, a beryllium window, filters, and a detector. The target angle was set to 20 
degrees. Monochromatic electrons were emitted from the electron gun and hit the 
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target to generate photons. The photons propagate through the beryllium window 
and the filtration material and arrive at the detector. The filtration material used in 
this work was aluminium. Photon detection was achieved by recording photons (both 
bremsstrahlung and characteristic X-rays) using the detector located at 300 mm from the 
anode. The size of the detector was limited to avoid oblique photon events. Aluminium 
filters of different thicknesses were added to simulate a set of model spectra. 

A sensitive detector was used to store information regarding interactions of a 
particle, and was attached to the detector volume to extract and store photon energy. 
A so-called “sensitive detector” is a detector model in GEANT4 which tracks particle 
interactions. Instead of recording the deposited energy of photons within the detector, 
we directly counted photons and extracted their energy using the sensitive detector. 
Thus the Geant4-based MC simulation can provide an accurate energy spectrum of the 
x-ray source by excluding the contribution of the detector, such as K-edge effect of the 
detector material. For each of the x-ray model spectrum, a total of 5 x 10 8 electrons were 
simulated. ROOT, an object oriented data analysis framework was used to analysis the 
recorded photon events. Energy spectrum of x-ray photons for different thicknesses of 
filtration were recorded as histograms. 

In order to estimate the spectrum as accurately as possible, it is recommended that 
the model spectra should cover the unknown true spectrum. Specifically, the energy 
range of the model spectra should be at least as broad as that of the true spectrum. For 
example, if the true spectrum was a 70 kVp polychromatic spectrum filtered with 5 mm 
of aluminium, then the softest model spectrum should be softer than the true spectrum 
and the hardest model spectrum should be harder than the true spectrum. 

Even though MC simulation can generate unattenuated model spectra which 
guarantee the actual raw spectrum will be covered, it is not a necessary requirement 
to generate sufficient model spectra with MC simulations. Other spectrum generators, 


such as SpekCalc (Poludniowski and Evans, 2007; Poludniowski, 2007; Poludniowski 


et a/., 2009), can also be used to generate a set of model spectra with a large dynamic 


range. 


3. Implementation 

The proposed spectrum estimation method was evaluated using both simulated and 
experimental data. In the following subsections, we will present the implementation 
details. 


3.1. Simulation study 

3.1.1. Simulation settings In order to make the numerical simulations as realistic 
as possible, a numerical phantom was created using a phantom image acquired on a 
cone beam CT scanner. The acquired water tank image (figure [2j) was first segmented 
to identify the water filling and the PMMA boundary. The patient table was also 















Figure 2. (a) Numerical phantom, (b) the 80 kVp and 100 kVp incident spectra 
generated using Spektr software. 


recognized as PMMA. Raw projection data was generated by polychromatic forward 
projecting the numerical phantom with an analytical calculated spectrum, purposely 
differing from the model spectra that will be used in the calibration procedure. Two 
sets of analytical polychromatic spectra with tube potentials of 100 and 80 kVp were 


generated using the spectrum generator Spektr (Siewerdsen et al ., 2004) and both 
of them were filtered with 3.5 mm aluminium. Since the raw projection data were 
generated using these two analytical spectra, the two spectra can serve as ground truth 
for evaluation. 

The source to isocenter distance of the simulated system was 1000 mm and the 
source to detector distance was 1500 mm. A circular scan which contains 675 projections 
in 360 degree angular range was simulated. Only one row of detector pixels (the central 
row) was simulated. The pixel size was 0.776 mm and the detector had 504 pixels. The 
propagation path length (PPL) for each material for each pixel was calculated using a 
GPU-based ray tracing algorithm, as illustrated in figure [3} Projections for each pixel 
were then calculated with the 100 kVp and 80 kVp analytical spectra and the PPL 
using the polychromatic projection formula (equation (|2j)) . After taking the logarithmic 
operation, the projection was regarded as raw projection data. During the spectrum 
estimation, all of the model spectra are normalized and the softest and the hardest 
model spectra always cover the raw spectrum for all evaluations. Because we can not 
yield the correct spectrum if the raw spectrum is out the dynamic range of the model 
spectra, regardless of how we combine the model spectra. During the weight calibration 
procedure, all weights were initialized with the same value. 


3.1.2. Number of model spectra In our first evaluation, the number of model spectra 
that were needed to yield an accurate estimated spectrum was determined and the 
influence of the number of model spectra on the accuracy of the estimated spectrum 
was also investigated. For the 100 and 80 kVp raw spectra, spectrum estimations with 
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Figure 3. Propagation path length in PMMA (a) and water (b) of the numerical 
water tank phantom. 

2, 4, and 6 model spectra were performed. 

In order to quantitatively characterize the accuracy of the estimated spectrum, we 
use the normalized root mean square error NRMSE and the mean energy difference 
between the true spectrum and the estimated spectrum A E which were defined as 
follows: 


NRMSE = 


Ef =1 (0(e)-^(e)) 2 

ZlMe ) 2 


( 10 ) 


N 

AE = J2 E e(.tt{e) -tl(e)) 

e=l 


( 11 ) 


Here Q(e) is the eth energy bin of the normalized estimated spectrum and Q(e) is eth 
energy bin of the normalized raw spectrum. N is the total number of the energy bins 
of the spectrum, and E(e) is the energy of the eth energy bin. 


3.1.3. Number of phantom materials The aim of the proposed method is to estimate 
the raw spectrum using the corresponding raw projection data. Nevertheless, in realistic 
applications, it is necessary to handle projections of objects with a varying number of 
materials. Thus, to make the method as realistic as possible, we also investigated the 
impact of the number of phantom materials on the accuracy of the estimated spectrum. 
The numerical phantom was generated using CT images of a water tank. In an initial 
simulation, both the water filling and the PMMA boundary were considered to be water, 
thus the numerical phantom only has one material. Then the raw projection data 
was generated by projecting the numerical phantom with the polychromatic projection 
formula. In a second simulation, the PMMA boundary was maintained and the raw 
projection data was generated by a polychromatic projection of the two materials. The 
raw spectra were estimated with six model spectra in both of the simulations. The 
accuracy of the estimated spectra was also evaluated with equations (10) and 0 - 
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3 .I. 4 . Robustness to model spectra generator The model spectra used in the previous 
simulations were generated using a MC tool-kit. However, compared to other spectrum 
generators, this process was time consuming and a large number of photon events are 
needed to yield a spectrum with good statistical properties. Other spectrum calculators, 


such as SpekCalc (Poludniowski et al, 2009), are much esier to handle and are well 


validated. Spectrum generation using these generators is very fast (usually on the order 
of seconds) and this time savings may play an important role in some applications, such 


as real-time full-spectral image reconstruction (Cai et al., 2013). Therefore, we also 


investigated the robustness of the proposed method with respect to the model spectra 
generator. 

For the numerical water tank phantom, a 70 kVp analytical raw spectrum which 


was modelled using Spektr (Siewerdsen et al, 2004), was employed to yield the raw 


polychromatic projection data. In the spectrum estimation procedure, the model 
spectra which were generated using both Geant4 and SpekCalc were used to estimate 
the raw spectrum. The number of model spectra were selected according to the 
previous simulation results. The accuracy of the estimated spectra were quantitatively 
characterized using NRMSE and A E. 

3.2. Experimental study 

The proposed method was also evaluated using experimental phantom data. The image 
uniformity module of the Catphan® G0(J/] phantom which was cast from a uniform 
material, was used to estimate the corresponding spectra of the acquisition protocols of 
the CT scanner. The raw projection data was acquired with a clinical bi-plane C-Arm 
angiography system (Artis zee, Siemens AG, Forchheim, Germany). Table [l] shows an 
overview of the system parameters and acquisition parameters. 

The C-Arm scanners were equipped with auto exposure control (AEC) systems 
which can perform automatic angular modulation of the tube current, tube potential or 
both. The modulation was performed according to the object’s size, shape and materials 
distribution to enhance radiation dose efficiency. For the Catphan®600 phantom CT 
scan, the tube potential was modulated between 82.5 kV and 90 kV. In this case, we 
can either estimate the spectrum of each view angle by using the corresponding raw 
projection data of each view angle, or alternatively, estimate an effective spectrum 
for the whole scan by using raw projections of all view angles. The accuracy of the 


effective spectrum can be evaluated using full spectrum reconstruction (Cai et al, 2013) 
or artefacts correction to see whether there are residual artefacts. However, this is 
out the scope of this study. In this study, the spectrum of a single view angle (the 
corresponding tube potential is 83.8 kV) was estimated independently, and its accuracy 
was evaluated by comparing the estimated spectrum to the spectrum generated using 


SpekCalc software (Poludniowski et al.\ 2009) with a beam hardening equivalent of 2.5 


mm of aluminium filtration (matching the value provided by the scanner manufacturer). 


t http://www.phantomlab.coin/library/pdf/catphan600_download.pdf. 
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Table 1 . System parameters and acquisition parameters of the experimental study. 


Parameters 

Artis zee 

Source to detector distance 

1200 mm 

Source to isocenter distance 

749.3 mm 

Angular range 

198.1 

View-angle increment 

0.4 

Number of views 

496 

X-ray source 

82.5 kVp - 90 kVp 

Cone angle 

00 

o 

Scanning time 

20 s 

Detector size 

380 x 296 mm 2 

Number of detector pixels 

1232 x 960 (with 2x2 rebinning) 

Detector material 

aSi with Csl scintillator 


To reduce the computational complexity with no loss of generality, only the central 
row of the raw projection data was extracted and compared to the polychromatic 
reprojection data which was calculated with a projection matrix. The raw spectrum 
was recovered by updating the weights, c, to minimize the quadratic error between the 
raw projection data and the reprojection data. During the weights calibration procedure, 
all of the initial guess were set to the hardest spectrum of the model spectra and the 
standard NIST attenuation coefficients of water was assigned to the segmented uniform 
module. 

Since the logarithmic raw projection data will be compared to the polychromatic 
reprojection data which doesn’t model the scatter signal, the scatter radiation 
contribution should be minimized in the raw projection data. In the experimental 
evaluation, the scanner was equipped with a focused anti-scatter grid and a narrow 
collimation mode was employed during the CT scan. Therefore, the contribution of the 
scatter signal to the raw projection was neglected. 

Compared to the diagnostic CT scanner, the absorption efficiency of the flat 
detector equipped on the C-Arm scanner is relatively low, especially for high energy 
photons. In this evaluation, a spectrum estimation both with and without taking the 
absorption efficiency into account was performed. After incorporating the absorption 
efficiency, the energy-dependent response of an energy-integrating detector should be 
the multiplication of the photon energy and the absorption efficiency. The energy 
dependent absorption efficiency was calculated using both the analytical Beer’s law 
and MC simulations with different energy levels. 
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(a) (b) 



Figure 4. 80 (a) and 100 (b) kVp model spectra generated using Geant4 with different 
thickness of aluminium filtration. 

Table 2. Weights of the model spectra of spectrum estimation for the numerical water 
tank phantom using six model spectra. 


Tube potential (kV) 



Weights 



Cl 

c 2 

c 3 c 4 

^5 

Cg 


80 

0.0730 

0.4023 

0.2874 0.1561 

0.0493 

0.0319 

100 

0.0171 

0.0369 

0.9379 0.0038 

0.0023 

0.0020 


4. Results 

4-1. Numerical phantom 

4 -1.1. Number of model spectra Figure [4] shows the 80 kVp and 100 kVp model spectra 
that were used in the numerical evaluation. These model spectra were generated using 
MC simulation tool-kit Geant4 with different thicknesses of aluminium filters. Figure [5] 
shows the results of the estimated 80 kVp and 100 kVp spectra as a function of the 
number of MC model spectra, using the numerical water tank phantom. All mean 
energy differences, A E, and normalized root mean square errors, NRMSE, are derived by 
comparing the estimated spectra with the raw analytical spectra depicted in figure [2j^b) . 

Figure [5^al), (a2) and (a3) show the 80 kVp spectra estimated with 2, 4 and 6 
model spectra and figure [5](bl) , (b2) and (b3) show the 100 kVp spectra estimated with 
2, 4 and 6 model spectra, respectively. For both of the two cases, NRMSE decreased as 
the number of model spectra increased, while A E didn’t have a significant change with 
a different number of model spectra. As can be seen, the spectrum estimation with 6 
model spectra leads to minor improvements in the accuracy. However, estimation with 
four model spectra can still yield acceptable results. Table [2] shows the calibrated weights 
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(al) 



(a2) 



(a3) 



(bl) 


-raw 


ZlE=0.25keV 

- estimated 


NRMSE=16% 



Energy (keV) 


02 ) 
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Figure 5. X-ray spectra estimated from the numerical phantom using different number 
of model spectra, (al), (a2) and (a3) are the 80 kVp spectra estimated using 2, 4 and 6 
Monte Carlo model spectra, respectively, (bl), (b2) and (b3) are the 100 kVp spectra 
estimated using 2, 4 and 6 Monte Carlo model spectra, respectively. The difference 
between the estimated spectra and the raw spectra are marginal when six model spectra 
are used. 
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of the model spectra of spectrum estimation for the numerical water tank phantom using 
six model spectra. 

NRMSE of the 100 kVp estimated spectra are much larger than NRMSE of 
the 80 kVp estimated spectra, even though both estimated spectra match the raw 
analytical spectra quite well. This discrepancy can be attributed to the mismatch of 
the characteristic peak which is of larger magnitude for the 100 kVp spectrum than 
the 80 kVp spectrum. The model spectra were generated using Geant4-based MC 
simulation tool-kit and were binned to 1 keV intervals. However, the raw spectra were 


generated using the Spektr software which is based on the TASMIP model (Boone and 


Seibert, 

1997; 

Sisniega et al, 

2013 

)and the TASMIP model is further based completely 


on interpolation of Fewell’s direct spectroscopy measurements results (Fewell et al. 


1981), which were tabulated at intervals of 2 keV. Thus, even though the raw Spektr 


values were interpolated to 1 keV intervals, it is difficult for the estimated spectrum to 
match the raw analytical spectra at the characteristic peak exactly. 

Based on these results, we chose six model spectra such that the NRMSE was 
minimized for all further simulations. It has to note that the spectra shown in figure [4] 
and [5] are estimated using the sequential optimization approach, however, both BFGS-B 
and ALM methods can provide the similar results. In this study, we use the sequential 
optimization approach for the sake of simplicity. 


4-1.2. Number of phantom materials In figure [6j the influence of the number of 
phantom materials used in estimation procedure on accuracy of the estimated spectrum 
with the proposed method is shown. Figure |6](a) and (b) depict the estimated 80 kVp 
and 100 kVp spectra with only one material. For better demonstration, both estimations 
use six MC model spectra. Compared to figure [5](a3) and (b3) which also use six model 
spectra to estimate the 80 kVp and 100 kVp spectra, the NRMSE increase from 4% to 6% 
and from 11% to 13% for the 80 kVp and the 100 kVp spectra, respectively. The reason 
for this behaviour can be attributed to the attenuation property of the polychromatic 
spectrum within different kinds of materials. If a single material phantom is used in 
the spectrum estimation, the estimated spectrum only takes the attenuation property 
of one material into account during the spectrum recovery. However, when a phantom 
consisting of more than one material is used, the estimated spectrum has to take all of 
the attenuation properties of these materials into account. This makes the spectrum 
recovery problem more deterministic. In practice, multi-materials phantom or a multi¬ 
materials CT image is recommended for the spectrum estimation problem. 


4-1.3. Robustness to model spectra generator In order to evaluate the robustness of the 
proposed method with respect to model spectra generator, both Geant4 model spectra 
and SpekCalc model spectra were used to recover the raw Spektr analytical spectrum. 
In figure [7|a), a set of 70 kVp model spectra which were generated using MC tool¬ 
kit Geant4 with different thickness of filtration were used in the weights calibration 
procedure. In figure [7](b) , the model spectra were generated using SpekCalc. The initial 
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(a) (b) 



Figure 6. Both 80 (a) and 100 (b) kVp x-ray spectra estimated from the numerical 
water tank phantom, assuming a single object material by setting the PMMA boundary 
to water. Compared to the spectra estimated using two materials, the spectra 
estimated using one material show inferior results. 

(a) (b) 




Figure 7. X-ray spectra estimated for the numerical water tank phantom data. 
Model spectra generated using both Geant4 (a) and SpekCalc (b) were employed in 
the estimation procedure to demonstrate the proposed method was robust against 
the model spectra generator. The initial guesses for the spectrum recovery problem 
correspond to the hardest model spectrum. Spectra estimated using both model 
spectra yield acceptable results, showing the proposed method was robust with respect 
to the model spectra generator. 


guess of the iterative calibration procedure was the hardest spectrum of the model 
spectra. Quantitative analysis of NRMSE and A E show spectra estimated using both 
model spectra generators can yield acceptable results, indicating the proposed method 
is robust with respect to the model spectra generator. Note that all of the spectra have 
1 keV energy bin widths and are normalized. 
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Figure 8. Energy dependent absorption efficiency of the Csl scintillator calculated 
using both the analytical method and MC simulation. Usually, the absorption efficiency 
goes down as photon energy increases, the absorption efficiency around 100% for low 
energies, however the K-edges of Cs and I elevate the absorption efficiency at 36.0 keV 
and 33.2 keV, respectively. 


4-2. Experimental phantom data 

Spectrum estimation using the proposed method with physical phantom was also 
performed on a clinical bi-plane C-Arm system (Artis zee, Siemens AG, Forchheim, 
Germany). Because of the low absorption efficiency of the flat detector that was 
equipped on the C-Arm scanner, we need to incorporate the absorption efficiency in the 
polychromatic reprojection procedure. Figure [8] shows the energy dependent absorption 
efficiency of the flat detector with 0.6 mm thickness of Csl. Both analytical calculation 
and MC simulation results are depicted. Usually, the absorption efficiency goes down 
as photon energy increases. However, the K-edges of Cs and I elevates the absorption 
efficiency at 36.0 keV and 33.2 keV, respectively. 

By taking the energy dependent absorption efficiency into account, the detector 
response is described as the multiplication of the photon energy and the absorption 
efficiency for an energy-integrating detector. Figure [9] shows the estimated 83.8 kVp 
spectrum by using the Catphan@600 phantom data with and without the incorporation 
of absorption efficiency. For this experimental phantom study, all of the model spectra 
used in the weights calibration, were generated by the SpekCalc software. Figure |9](a) 
shows the spectrum estimated without taking the absorption efficiency into account. 
The reference raw spectrum was generated using SpekCalc software with x-ray tube 
specifications provided by the scanner manufacturer. The initial guess was the hardest 
model spectrum. As can be seen, the mean energy of the estimated spectrum was 1.25 
keV lower than the raw spectrum. This was because the experimental raw projection 
data tended to yield a lower photon count due to the low absorption efficiency. However, 
the absorption efficiency was not taken into account (i.e., it was regarded as 100%) in the 
polychromatic reprojection procedure. Thus, to compensate the difference of the raw 
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(a) (b) 



Figure 9. X-ray spectra estimated for the C-Arm system using the Catphan®600 
phantom data with (a) and without (b) absorption efficiency incorporation. Model 
spectra generated using SpekCalc were employed in the estimation procedure, and the 
raw spectrum was generated using SpekCalc software and x-ray tube specifications 
provided by the scanner manufacturer. The initial guesses of the spectrum recovery 
problem correspond to the hardest model spectrum. Due to the low absorption 
efficiency, the estimated spectrum can not match the raw spectrum well without taking 
the absorption efficiency into account for flat panel detector based cone-beam CT, 
but by incorporating the absorption efficiency, the estimated spectrum can be greatly 
improved. 


projection data and the polychromatic reprojection data, the estimated spectrum tended 
to be softer. By taking the photon absorption efficiency into account (figure [9](b)), the 
estimated spectrum matched the raw spectrum much better. Both of NRMSE and A E 
were greatly reduced. 

5. Discussion 

The proposed method was based on transmission measurement, namely, by updating the 
estimated spectrum to minimize the difference between the raw projection data and the 
polychromatic reprojection data. In this study, the estimated spectrum was expressed 
as a weighted summation of a set of model spectra. Instead of directly estimating the 
content of each energy bin of the spectrum, the proposed spectrum estimation method 
estimates the corresponding weights for each model spectrum. Thus the number of 
unknown variables for the spectrum recovery problem can be greatly reduced from the 
number of energy bins to the number of model spectra and the solution of inverse 
problem was very stable, especially when two materials phantom were employed in 
the weights calibration procedure. Quantitative evaluations show the proposed method 
yields acceptable spectra with four model spectra, and the accuracy may be further 
improved with an increased number of model spectra. Since the proposed method has 
shown robustness with respect to model spectra generator and it is very convenient to 
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generate spectra with different filtrations using SpekCalc, it is recommended to use a 
higher number of model spectra in the spectrum recovery procedure, such as six or eight 
model spectra. 

The rationale of a higher number of model spectra yield improved estimation results, 
is that the model spectra make up a database and the spectrum calibration procedure 
is to find an optimal spectrum using the database to generate a reprojection dataset 
whose difference is as small as possible from the raw projection data. Thus the accuracy 
of estimated spectrum could improve as the sufficiency of the database increases. When 
only two model spectra(the softest and the hardest model spectra, i.e., a small database) 
were used, it is difficult for these two spectra to make up a medium filtered spectrum if 
the raw spectrum is a medium filtered spectrum, as shown in figure [5](al) and (bl). In 
this case, the estimated spectrum is inferior due to the insufficiency of the model spectra 
database. When six model spectra were used, the estimated spectra were significantly 
improved since it is much easier to find a spectrum using the database which has a 
similar medium filtered profile (shown in figure [5](a3) and (b3)). One may want to use 
a large amount of model spectra (a large database) in the weights calibration, however, 
the calibration time may increase as the number of model spectra increases. 

In order to keep the model spectra covering the incident spectrum, one may use 
a broad filtration during the model spectra generation. It has to note that since the 
spectrum is the weighted summation of a set of model spectra, the model spectra do 


impact on the final result. For example, spectrum estimation using Spektr (Siewerdsen 


et al ., 2004) model spectra and SpekCalc (Poludniowski et a/., 2009) model spectra may 


have different profiles, even though both of the results may have similar mean energies 
and have minimal NRMSEs. Because most of the widely used spectrum generators are 
well verified, spectrum estimation using these generators would yield acceptable results. 

Since the proposed method tries to minimize the quadratic error of the raw 
projection data and the estimated reprojection data. The contribution of the raw 
projection data includes not only the attenuation of the object, but also the extra 
attenuation (such as bow-tie filter, detector sensor protection materials). Thus the 
estimated spectrum which was used to calculate the estimated reprojection data, should 
be an effective spectrum which has taken these extra attenuation into account. In 
addition, the method described here is localized; in other words, the estimated spectrum 
corresponds to those detector pixels which provide its raw projection data for the 
spectrum estimate. Since the central detector row of the raw projection data was used 
to calibrate the spectrum, the estimated spectrum can be considered to be an effective 
spectrum for the whole fan beam. The heel effect of the x-ray tube is not taken into 
account here. However, to estimate an average spectrum along different cone angle 
which incorporates the average attenuation of the heel effect, raw projection data of the 
central column of the detector may be employed. Although the spectra along different 
fan angle and cone angle are different, an average spectrum or effective spectrum may 
be sufficient for most the realistic spectrum applications. 

Because of its pixel-wise feature, the proposed method can be applied to estimate 
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the spectrum of an CT scanner equipped with auto exposure control (AEC) system. To 
enhance radiation dose efficiency, CT scan with AEC system perform automatic angular 
modulation of the tube current, tube potential or both from view by view by according to 
the object’s size, shape and materials distribution. Thus the x-ray spectrum may change 
from view by view during the CT scan. To estimate the spectrum of a specific view 
angle, the raw projection data of this specific view angle can be employed to estimate 
the corresponding spectrum. In addition, by using the model spectra of different tube 
potential and projections from all of the view angles, the proposed method can be further 
applied to estimate an effective spectrum for a whole CT scan with tube potential 
modulation. In this case, the correctness of the estimated effective spectrum has to 
be evaluated using either a full-spectrum image reconstruction, artefact correction or 
weighted spectra for all view angles. Our aim here was not to investigate the correctness 
of the estimated effective spectrum but to evaluate the proposed spectrum estimation 
method. Therefore, the evaluation of the correctness of the effective spectrum of all view 
angles is not within the scope of this paper. In future studies, we will investigate how 
to estimate an effective spectrum and to evaluate its correctness using either artefacts 
correction or weighted spectra from all view angles. 

It is important to exclude scattered radiation in the raw projection data. The 
proposed process compares the raw projection data with the polychromatic reprojection 
data which does not contain scatter signal. The experimental data based on C-Arm 
measurements arguably contains a comparably high scatter fraction. However, this did 
not show considerable influence on the results. This might be attribute to the focused 
scatter grid employed by the systems, rejecting a considerable amount of scattered 
radiation. Further work is necessary to test the proposed method to the influence of 
scattered radiation. 

Compared to the diagnostic CT scanner, the energy dependent absorption efficiency 
of the fiat-panel detectors (FPD) that were used in the C-Arm systems is relatively low 
because the stopping power of Csl-based FPD pixel is inferior than the Gd 2 C> 2 S-based 
diagnostic CT detector pixel and the thickness of the Csl scintillator is also limited. 
This is the reason that there was clear discrepancy between the spectrum estimated with 
and without photon absorption efficiency incorporation for the experimental evaluation. 
Since the absorption efficiency is much higher for the diagnostic CT scanner, to our 
belief, even without taking the absorption efficiency into account, spectrum estimation 
for the diagnostic CT scanner may yield acceptable result. 

6. Conclusions 

We proposed an indirect transmission measurement-based spectrum estimation method 
that uses the raw projection data of a phantom with known density to recover the 
corresponding effective spectrum. The method has been evaluated using both simulated 
data and experimental data. Quantitative analysis demonstrated the method can 
estimate the x-ray spectra quite well. Further tests showed that a higher number of 
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model spectra and multi-component phantoms or objects can yield improved spectrum 
estimation. In addition, the proposed method is robust with respect to the model spectra 
generator, thus several widely used spectra generators (such as Spektr and SpekCalc) 
can be employed to estimate the spectrum. 
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